﻿// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

using System;
using System.Linq;
using Microsoft.Quantum.Simulation.Core;
using Microsoft.Quantum.Simulation.Simulators;
using Microsoft.Quantum.Simulation.Simulators.QCTraceSimulators;
using Microsoft.Quantum.Chemistry;
using Microsoft.Extensions.Logging;
using System.Diagnostics;
using Microsoft.Quantum.Chemistry.OrbitalIntegrals;
using Microsoft.Quantum.Chemistry.Fermion;
using Microsoft.Quantum.Chemistry.QSharpFormat;

// This loads a Hamiltonian from file and performs quantum phase estimation on
// - Jordan–Wigner Trotter step
// - Jordan–Wigner optimized Trotter step
// - Jordan–Wigner Qubitization iterate

namespace Microsoft.Quantum.Chemistry.Samples
{

    class Program
    {

        static void Main(string[] args)
        {
            var logger = Logging.LoggerFactory.CreateLogger<Program>();

            // Directory containing integral data generated by Microsoft.
            //Example Liquid data format files
            /*
            "h2_sto3g_4.dat" // 4 SO
            "B_sto6g.dat" // 10 SO
            "Be_sto6g_10.dat" // 10 SO
            "h2o_sto6g_14.dat" // 14 SO
            "h2s_sto6g_22.dat" // 22 SO
            "co2_sto3g_30.dat" // 30 SO
            "co2_p321_54.dat" // 54 SO
            "fe2s2_sto3g.dat" // 110 SO
            "nitrogenase_tzvp_54.dat" // 108 SO
            */
            


            // Read Hamiltonian terms from file.
            // Stopwatch for logging time to process file.
            Stopwatch stopWatch = new Stopwatch();


            #region For loading data in the format consumed by Liquid.
            stopWatch.Start();
            string LiquidRoot = @"..\IntegralData\Liquid\";
            string LiquidFilename = "h2_sto3g_4.dat";
            
            logger.LogInformation($"Processing {LiquidFilename}");
            var problemDescriptionLiQuiD = LiQuiD.Deserialize($@"{LiquidRoot}\{LiquidFilename}").Single();
            // Number of electrons. This must be specified for the liquid format.
            problemDescriptionLiQuiD.NElectrons = 2;

            logger.LogInformation($"Liquid Load took {stopWatch.ElapsedMilliseconds}ms");
            stopWatch.Restart();
            #endregion
            #region For loading data in the Broombridge format.
            stopWatch.Start();
            string YAMLRoot = @"..\IntegralData\YAML\";
            string YAMLFilename = "lih_sto-3g_0.800_int.yaml";
            var problemDescriptionBroombridge = Broombridge.Deserializers.DeserializeBroombridge($@"{YAMLRoot}\{YAMLFilename}").ProblemDescriptions.Single();
            logger.LogInformation($"Broombridge Load took {stopWatch.ElapsedMilliseconds}ms");
            stopWatch.Restart();
            #endregion

            // Select problem description to use
            var problemDescription = problemDescriptionBroombridge;
            var fermionHamiltonian = problemDescription.OrbitalIntegralHamiltonian.ToFermionHamiltonian(IndexConvention.UpDown);

            // Logs spin-orbital data in Logger.Message.
            logger.LogInformation(fermionHamiltonian.ToString());

            // Process Hamiltonitn to obtain Jordan–Wigner representation.
            var jordanWignerEncoding = fermionHamiltonian.ToPauliHamiltonian(Paulis.QubitEncoding.JordanWigner);

            // Create input wavefunction.
            var wavefunction = fermionHamiltonian.CreateHartreeFockState(problemDescription.NElectrons);

            // Alternately, use wavefunction contained in problem description,
            // if available
            // var wavefunction = problemDescription.Wavefunctions["label"].ToIndexing(IndexConvention.UpDown);

            // Logs Jordan–Wigner representation data in Logger.Message.
            logger.LogInformation(jordanWignerEncoding.ToString());

            logger.LogInformation("End read file");

            // We begin by making an instance of the simulator that we will use to run our Q# code.
            using (var qsim = new QuantumSimulator())
            {
                // Package hamiltonian and wavefunction data into a format
                // consumed by Q#.
                var qSharpData = QSharpFormat.Convert.ToQSharpFormat(
                    jordanWignerEncoding.ToQSharpFormat(),
                    wavefunction.ToQSharpFormat());

                #region Calling into Q#
                // Bits of precision in phase estimation.
                var bits = 10;

                // Repetitions to find minimum energy.
                var reps = 5;

                // Run phase estimation simulation using Jordan–Wigner Trotterization.
                var runJW = true;
                
                // Trotter step size.
                var trotterStep = .4;

                // Run phase estimation simulation using Jordan–Wigner Trotterization with optimzied circuit.
                var runJWOptimized = true;
                
                // Run phase estimation simulation using Jordan–Wigner qubitization.
                var runJWQubitization = true;

                if (runJW)
                {
                    for (int i = 0; i < reps; i++)
                    {
                        var (phaseEst, energyEst) = TrotterEstimateEnergy.Run(qsim, qSharpData, bits, trotterStep).Result;
                        logger.LogInformation($"Trotter simulation. phase: {phaseEst}; energy {energyEst}");
                    }
                }
                if (runJWOptimized)
                {
                    for (int i = 0; i < reps; i++)
                    {
                        var (phaseEst, energyEst) = OptimizedTrotterEstimateEnergy.Run(qsim, qSharpData, bits - 1, trotterStep).Result;
                        logger.LogInformation($"Optimized Trotter simulation. phase {phaseEst}; energy {energyEst}");
                    }
                }
                if (runJWQubitization)
                {
                    for (int i = 0; i < reps; i++)
                    {
                        var (phaseEst, energyEst) = QubitizationEstimateEnergy.Run(qsim, qSharpData, bits - 2).Result;
                        logger.LogInformation($"Qubitization simulation. phase: {phaseEst}; energy {energyEst}");
                    }
                }

                #endregion
            }


            if (System.Diagnostics.Debugger.IsAttached)
            {
                System.Console.ReadLine();
            }

        }
    }
}

